%% load data 

clear C S cell 

%passive and active together, passive = 1:10, active 11:40 (a1, a2, a30 
N = 10;
%C = cell(1,N);

bname = " "; %same file name as from CE_run

for k = 1:N
   path = " " + bname + "_" + num2str(k) + "/wspace.mat"; % input root folder that base name is in
   C{k} = load(path,'data');
   
end
S = [C{:}]; % create structure array


%% Gather data S(sim).data(datum, cell, time)


sims = [1:5; 5:10]; %number of sims, columns to compare groups


areas = [];
sideds = [];

for i=1:10

area = S(i).data(37,:,1);
areas = [areas, area];

sided = S(i).data(40,:,1);
sideds = [sideds, sided];

end

ep = [2:7, 8:15, 19]; %cell and corona endpt strains, PCD, MSD 

% find last time step with data 
hightime = 0; 
for j =1:numel(sims)
    timez = size(S(sims(j)).data,3);
    if (timez > hightime)
        hightime = timez; 
    end
end

% gather endpoint data 
endpt = {}; 

for j = 1:size(sims,1) %groups of sims (passive, active)
    for k=ep %across data types (strains, msd) 
    things = [];
        for i=sims(j,:)
            
                thing = S(i).data(k,:,:);
                thing = squeeze(thing);
             
                if (size(thing,2)) < hightime
                    thing = [thing, zeros(length(thing),hightime-size(thing,2))];
                end
                
        things = [things; thing];
        things(things==0)=NaN;
        end     
        
        endpt(find(ep==k),j)= {things};    
    end
end

 
 col =3;
 
 % Define 8 vectors of different lengths
v1 = ends{1,col};
v2 = ends{2,col};
v3 = ends{3,col};
v4 = ends{4,col};
v5 = ends{5,col};
v6 = ends{6,col};
v7 = ends{7,col};
v8 = ends{8,col};

% Store vectors in a cell array
vectors = {v1, v2, v3, v4, v5, v6, v7, v8};

% Find the maximum length of the vectors
max_length = max(cellfun(@length, vectors));

% Create a padded matrix with NaNs
padded_matrix = nan(length(vectors), max_length);

% Populate the matrix with the vectors
for i = 1:length(vectors)
    vector = vectors{i};
    padded_matrix(i, 1:length(vector)) = vector;
end
probs = padded_matrix';


 %getting timestep of stoppage from ends
 
  tsend = {};
  
  for i=1:size(sims,1)
      frz = []; 
      for j=1:size(sims,2)
          stop = size(S(sims(i,j)).data,3);
          frz = [stop; frz]; 
      end
      tsend{1,i} = frz; 
  end
  
  x = tsend{1,7}; 
  dx = [1, diff(x)];
x2 = x(dx ~= 0);

 tsend = {};
 
 for i=1:length(ends)
    tz = ends{1,i};
    tz = (tz-200)/3;
    tsend{1,i} = tz;
 end
 
 
 ends = {};
 for i=[1:15]  %just strains 1:6, 15 is msd, 7:14 are tx probabilities
     vals = [];
     for j=1:size(endpt,2)
         A = endpt{i,j};
         A(~any(~isnan(A), 2),:)=[]; %fix for txn probs with all NAN rows
         B = ~isnan(A);
         Indices = arrayfun(@(x) find(B(x, :), 1, 'last'), 1:size(A, 1),'UniformOutput',false);
         Indices = cell2mat(Indices);
         for k=1:length(A)
             vals(k) = A(k, Indices(k));
         end
         %Probabilities bc diff numbers of nan row sremoved per sim
         ends{i,j} = vals; %makes i*j, want all in one, do for
         %ends{1,j}(i,:) = vals;
     end
 end
 
 %same but for keeping NaNs 
 endN = {};
 

 for i=7:14%just strains 1:6, 15 is msd, 7:14 are tx probabilities
     vals = [];
     for j=1:size(endpt,2)
         A = endpt{i,j};
        
         B = ~isnan(A);
         IndicesMatrix = NaN(size(A, 1), 1);
         
         for m = 1:size(A, 1)
             idx = find(B(m, :), 1, 'last');
             if ~isempty(idx)
                 IndicesMatrix(m) = idx;
             end
         end

         vals = NaN(size(A, 1), 1);
         for k = 1:size(A, 1)
             if ~isnan(IndicesMatrix(k))
                 vals(k) = A(k, IndicesMatrix(k));
             end
         end
                
         endN{1,j}(i,:) = vals;
     end
 end
 
 

% time series data for step strains 
ts = [29:34, 20]; %inginal time series for concordance and t1s
ts = [37,40]; %area and sidedness sidedness is 40, area is 37
%fields

hightime = 0; 
for j =1:numel(sims)
    timez = size(S(sims(j)).data,3);
    if (timez > hightime)
        hightime = timez; 
    end
end

times = {}; 
for j = 1:size(sims,1) %groups of sims (passive, active)
    for k=ts %across data types (strains, msd) 
    things = [];
        for i=sims(j,:)
            thing = S(i).data(k,:,:);
            sthing = squeeze(thing);
            
            %fills in zeros for ones that didn't make it 
            if (size(sthing,2)) < hightime
                sthing = [sthing, zeros(length(sthing),hightime-size(sthing,2))];
            end
            
            things = [things;sthing]; %vertcat issue for diff times of stoppage, make sure hightime is used  
            things(things==0)=NaN; %off if want t1s w/o zeros for frequency
           
        end
        times(find(ts==k),j)= {things};    
    end
end

%% Further analysis of gathered data

%average of everything in endpoints
endavg = [];
for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(ends,1)
        endavg(m,j) = nanmean(ends{m,j});
    end
end

%get avg of end avg for analysis of weights
%includes ep strains, tx probs, msd 
end_all= endavg;


% do same for t1 dwell, freq, and leading using times 

% t1s 
t1_freq = nanmean(t1savg{1,1});

t1_dwell = nanmean(t1countavg{1,1}); 

%for Area +, just need area +
area_plus =mean(cfnsposavg{3,1}); 
 

% create a table of means, stds to calculate the metric for each case 
% metric variable in fig 7 folder, non zero values are z scores 
% rows:
% 1-6:Cell ML, Cell AP, Cor ML, Cor AP, Cell Area, Cor Area
% 7-8: T1 freq, dwell 
% 9: MSD
% 10-18: ml-, ml, ML+, AP-, ap, AP+, AREA-, AREA, AREA+
% 19-26: PDD ML, PDC ML, PCD ML, PCC ML, PDD AP, PDC AP, PCD AP, PCC AP

%columns 
% 1: P relative to AP, to calculate Ap score 
% 2: Ap relative to P, to calculate P score
% 3-5: contrac, relative to Ap, crawl, cap
% 6-8: Crawl, relative to Ap, crawl, cap
% 9-11: cap, relative to Ap, contrac, crawl

%calculate passive score, use column 1 
% sum of mean* weight over sum of weights

% make column of avgs that's the same length as metric
scores = [];
scores(1:6) = end_all(1:6);  
scores(7) = t1_freq;
scores(8) = t1_dwell;
scores(9) = end_all(15); 
scores(10:17) = 0;
scores(18) = area_plus;
scores(19:26) = end_all(7:14);

all_s = [];
all_s(1,:) =  ends{1,1}; 
all_s(2,:) =  ends{2,1};
all_s(3,:) = ends{3,1};
all_s(4,:) = ends{4,1};
all_s(5,:) = ends{5,1};
all_s(6,:) = ends{6,1};
all_s(7,:) = t1savgN';
all_s(8,:) = t1countavg{1,1};
all_s(9,:) =  ends{15,1}; 
all_s(10:17,:) = 0;
all_s(18,:) = cfnsposavg{3,1}';
all_s(19,:) = endN{1,1}(7,:); 
all_s(20,:) = endN{1,1}(8,:);
all_s(21,:) = endN{1,1}(9,:);
all_s(22,:) = endN{1,1}(10,:);
all_s(23,:) = endN{1,1}(11,:);
all_s(24,:) = endN{1,1}(12,:);
all_s(25,:) = endN{1,1}(13,:);
all_s(26,:) = endN{1,1}(14,:);


load('\metrics.mat') % add path
metrics = abs(metrics); % metrics are just weights, take abs
load('\means_ref.mat') % add path
load('\std_ref.mat') % add path

%means_ref is all means, std_refs is all stds 
%same rows as metrics, but cols are:
% 1. Passive, 2. Active, 3. contrac, 4. crawl, cap 


p_count = []; 

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,2))/std_ref(i,2); %z score of that item relative to the mean and std of reference, 
    %in this case active pooled, second col of mean and std ref matrices

    if metrics(i,1)~=0
        p_now =(item-metrics(i,1))^2;
        p_count = [p_count; p_now]; 
    end
end

p_score = sqrt(sum(p_count)) 

%for all cells
p_scores = [];

for i=1:size(all_s,2) %over cells 
    p_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,2))/std_ref(j,2); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,1)~=0
            now =(item-metrics(j,1))^2;
            p_counts = [p_counts; now];
        end
    end
    p_scores(i) = sqrt(nansum(p_counts, 'all'));
end

p_mean = mean(p_scores)


% active pooled similarity score, use passive z values as reference for
% distance 

ap_count = []; %[]; 
for i=1:length(metrics)
    item(i) = (scores(i)-means_ref(i,1))/std_ref(i,1); %z score of that item relative to the mean and std of reference, 

    if metrics(i,2)~=0  %only use sig metrics 
        p_now =(item(i)-metrics(i,2))^2;
        ap_count = [ap_count; p_now]; 
    end
end

ap_score = sqrt(sum(ap_count)) 

%ap score for all cells

ap_scores = [];

for i=1:size(all_s,2) %over cells 
    ap_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,1))/std_ref(j,1); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,2)~=0
            now =(item-metrics(j,2))^2;
            ap_counts = [ap_counts; now];
        end
    end
    ap_scores(i) = sqrt(nansum(ap_counts, 'all'));
end

ap_mean = mean(ap_scores)

% contraction, a1 score 

a1_count = []; %[]; 
% 3 refs for each active case

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,2))/std_ref(i,2);
    
    if metrics(i,3)~=0
        p_now =(item-metrics(i,3))^2;
        a1_count = [a1_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,4))/std_ref(i,4);
    
    if metrics(i,4)~=0
        p_now =(item-metrics(i,4))^2;
        a1_count = [a1_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,5))/std_ref(i,5);
    p_now =(item-metrics(i,5))^2;
    if metrics(i,5)~=0
        p_now =(item-metrics(i,5))^2;
        a1_count = [a1_count; p_now];
    end
end
    

a1_score = sqrt(sum(a1_count)) 

%all cells 
a1_scores1 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,2))/std_ref(j,2); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,3)~=0
            now =(item-metrics(j,3))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a1_scores1 = [a1_scores1; sqrt(nansum(a1_counts, 'all'))];
end

a1_scores2 = [];
for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,4))/std_ref(j,4); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,4)~=0
            now =(item-metrics(j,4))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a1_scores2 = [a1_scores2; sqrt(nansum(a1_counts, 'all'))];
end

a1_scores3 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,5))/std_ref(j,5); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,5)~=0
            now =(item-metrics(j,5))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a1_scores3 = [a1_scores3; sqrt(nansum(a1_counts, 'all'))];
end

%collect all a1 scores, one column per comp
a1_scores(:,1)= a1_scores1; % to Ap
a1_scores(:,2) = a1_scores2; % to crawl
a1_scores(:,3) = a1_scores3; % to cap

%means across rows
a1_means = mean(a1_scores'); 

% total mean
a1_mean = mean(a1_means)


% contraction a2 score 
a2_count = []; %[]; 
% 3 refs for each active case

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,2))/std_ref(i,2);
    
    if metrics(i,6)~=0
        p_now =(item-metrics(i,6))^2;
        a2_count = [a2_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,3))/std_ref(i,3);
    
    if metrics(i,7)~=0
        p_now =(item-metrics(i,7))^2;
        a2_count = [a2_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,5))/std_ref(i,5);
    if metrics(i,8)~=0
        p_now =(item-metrics(i,8))^2;
        a2_count = [a2_count; p_now];
    end
end
    

a2_score = sqrt(sum(a2_count)) 

%all cells crawl 
a2_scores1 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,2))/std_ref(j,2); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,6)~=0
            now =(item-metrics(j,6))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a2_scores1 = [a2_scores1; sqrt(nansum(a1_counts, 'all'))];
end

a2_scores2 = [];
for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,3))/std_ref(j,3); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,6)~=0
            now =(item-metrics(j,6))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a2_scores2 = [a2_scores2; sqrt(nansum(a1_counts, 'all'))];
end

a2_scores3 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,5))/std_ref(j,5); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,8)~=0
            now =(item-metrics(j,8))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a2_scores3 = [a2_scores3; sqrt(nansum(a1_counts, 'all'))];
end

%collect all a1 scores, one column per comp
a2_scores(:,1)= a2_scores1; % to Ap
a2_scores(:,2) = a2_scores2; % to crawl
a2_scores(:,3) = a2_scores3; % to cap

%means across rows
a2_means = mean(a2_scores'); 

% total mean
a2_mean = mean(a2_means)



% capture a3 score 
a3_count = []; %[]; 
% 3 refs for each active case

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,2))/std_ref(i,2);
    
    if metrics(i,9)~=0
        p_now =(item-metrics(i,9))^2;
        a3_count = [a3_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,3))/std_ref(i,3);
    
    if metrics(i,10)~=0
        p_now =(item-metrics(i,10))^2;
        a3_count = [a3_count; p_now];
    end
end

for i=1:length(metrics)
    item = (scores(i)-means_ref(i,4))/std_ref(i,4);
    if metrics(i,11)~=0
        p_now =(item-metrics(i,11))^2;
        a3_count = [a3_count; p_now];
    end
end
    

a3_score = sqrt(sum(a3_count)) 

%all cells contraction 
a3_scores1 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,2))/std_ref(j,2); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,9)~=0
            now =(item-metrics(j,9))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a3_scores1 = [a3_scores1; sqrt(nansum(a1_counts, 'all'))];
end

a3_scores2 = [];
for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,3))/std_ref(j,3); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,10)~=0
            now =(item-metrics(j,10))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a3_scores2 = [a3_scores2; sqrt(nansum(a1_counts, 'all'))];
end

a3_scores3 = [];

for i=1:size(all_s,2) %over cells 
    a1_counts = [];
    for j=1:length(metrics)
        item = (all_s(j,i)-means_ref(j,4))/std_ref(j,4); %z score of that item relative to the mean and std of reference,
        
        if metrics(j,11)~=0
            now =(item-metrics(j,11))^2;
            a1_counts = [a1_counts; now];
        end
    end
    a3_scores3 = [a3_scores3; sqrt(nansum(a1_counts, 'all'))];
end

%collect all a1 scores, one column per comp
a3_scores(:,1)= a3_scores1; % to Ap
a3_scores(:,2) = a3_scores2; % to crawl
a3_scores(:,3) = a3_scores3; % to cap

%means across rows
a3_means = mean(a3_scores'); 

% total mean
a3_mean = mean(a3_means)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LEADER ANALYSIS FOR STRAIN

% shifted corona strains

shifts=-5:5;
ms = [3,4,6]; %rows from step matrix with corona strains
shift = {};

for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:length(ms) %stored cor strain in {times} x, y, area
        cellshif = [];
        for i=sims(j,:)
            %get number of cells in this sim
            cells = size(S(i).data,2);
            cellshi = [];
            for k= 1:cells
                
                %get corona strains
                strain = times{ms(m),j}(k,:);
                cellsh = []; %shifts for one cell for one strain type
                for s=1:length(shifts)
                    col = circshift(strain, shifts(s));
                    if shifts(s)<0
                        col(shifts(s+1)+length(strain):length(strain)) = NaN;
                    elseif shifts(s)>0
                        col(1:shifts(s)) = NaN;
                    end
                    cellsh(s,:) = col; %shifts by time for one cell
                end
            cellshi(:,:,k) = cellsh;
            
            end
            cellshif = cat(3,cellshif,cellshi);
        end
        
        shift(m,j) = {cellshif};
    end
end

% metric

cs = [1,2,5]; %rows from step matrix with cell strains
cfns = {}; %dnt think this stores anything, see cnfs for that info

for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:length(cs) %stored cor strain in {times} x, y, area
        cfncelli = [];
        for i=sims(j,:)
            %get number of cells in this sim
            cells = size(S(i).data,2);
            cfncell = [];
            for k= 1:cells
                
                cellst = times{cs(m),j}(k,:);
                
                cellsh = [];
                for s=1:length(shifts)
                    
                    corst = shift{m,j}(s,:,k);
                    
                    
                    for t=1:length(cellst)
                        
                        cell = cellst(t);
                        cor = corst(t);
                        
                        
                        if ~isnan(cor) && ~isnan(cell)
                            
                            se = abs(cell - cor); % dist from avg-- same for either
                          
                            if (sign(cell) == sign(cor))
                                adj = (-1*se) +1;
                            else
                                adj = -1*se;
                             
                            end
                            
                            if (cell && cor)==0
                                adj=0;
                            end
                            
                           
                        else
                            adj= NaN;
                        end
                        
                        cellsh(s,t) = adj;    
                        
                    end
                    
                   
                end
            cfncell(:,:,k) = cellsh;     
            end
        cfncelli = cat(3, cfncelli, cfncell); 
        end
        
        cfns(m,j) = {cfncelli};
    end
end

%average over negatives
cfnsneg = {};


for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area

          cfnsneg{m,j} = cfns{m,j}(1:5,:,:);
    end

end

cfnsnegavg = {};

for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area
            avgs = nanmean(cfnsneg{m,j}(:,:,:), [1,2]);
            cfnsnegavg{m,j} = squeeze(avgs);
    end

end

% zeros
cfns0 = {};


for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area

          cfns0{m,j} = cfns{m,j}(5,:,:);
    end

end

cfns0avg = {};

for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area
            avgs = nanmean(cfns0{m,j}(:,:,:), [1,2]);
            cfns0avg{m,j} = squeeze(avgs);
    end

end


% average over positives
cfnspos = {};


for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area

          cfnspos{m,j} = cfns{m,j}(6:11,:,:);
    end

end


cfnsposavg = {};

for j = 1:size(sims,1) %groups of sims (passive, active)
    for m=1:size(cfns,1) %stored cor strain in {times} x, y, area
            avgs = nanmean(cfnspos{m,j}(:,:,:), [1,2]);
            cfnsposavg{m,j} = squeeze(avgs);
    end

end

call = {};

%neg first 
for i=1:size(sims,1) %column of sim groups, p, a1, a2, a3 
    for j=1:3 %ML, AP, area 
        it = cell2mat(cfnsnegavg(j,i))';
        k=[1,4,7];
        call{1,i}(k(j),:) = it; 
       
    end 
        
end

%0s
for i=1:size(sims,1) %column of sim groups, p, a1, a2, a3 
    for j=1:3 %ML, AP, area 
        k =[2,5,8];
        it = cell2mat(cfns0avg(j,i))';
        call{1,i}(k(j),:) = it;       
        
    end 
        
end

%pos
for i=1:size(sims,1) %column of sim groups, p, a1, a2, a3 
    for j=1:3 %ML, AP, area 
        k=[3,6,9];
        it = cell2mat(cfnsposavg(j,i))';
        call{1,i}(k(j),:) = it; 
        
    end 
        
end

%collect them into 9 columns (ML-, ML, ML+, AP-, AP, AP+, Area-, Area,
%Area+) per sim column

simcol = 3;
all9 = [];

all9(:,1) = cfnsnegavg{1,simcol}; %ML- 
all9(:,2) = cfns0avg{1,simcol};
all9(:,3) = cfnsposavg{1,simcol};

all9(:,4) = cfnsnegavg{2,simcol}; %AP- 
all9(:,5) = cfns0avg{2,simcol};
all9(:,6) = cfnsposavg{2,simcol};

all9(:,7) = cfnsnegavg{3,simcol}; %area-
all9(:,8) = cfns0avg{3,simcol};
all9(:,9) = cfnsposavg{3,simcol};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% t1 transition average 

t1savg = {}; %without NaNs

for j = 1:size(sims,1) %groups of sims (passive, active)
     %t1s are in row 7, each one is all cells x time, wanna average over time
     %get rid of nans
            use = times{7,j}; %or 7,j if all calculated 
          
            use(~any(~isnan(use), 2),:)=[];
            use  = ~isnan(use);
          
            avgs = nanmean(use'); 
            t1savg{1,j} = squeeze(avgs);
end

t1savgN = []; %with NaNs

for j = 1:size(sims,1) %groups of sims (passive, active)
    use = times{7,j}; %or 7,j if all calculated
    for i=1:size(times{7,1},1)
        tot = nansum(use(i,:),'all');  
        if tot == 0
            t1savgN(i,j) = NaN;
        else 
            t1savgN(i,j) = tot/size(times{7,1},2);
        end
    end
    
end



%t1 transition dwelling times


t1count = {}; 

for k = 1:size(sims,1) 
    
inputMatrix =  times{7,k}; %use 7 if step strains also in 'times' 
[numEntities, numTimepoints] = size(inputMatrix);


outputMatrix = zeros(numEntities, numTimepoints);


for i = 1:numEntities

    count = 0;

    colIndex = 1;
    
    % Loop through each timepoint 
    for j = 1:numTimepoints
        if inputMatrix(i, j) == 1
        
            count = count + 1;
        else
         
            if count > 0
                outputMatrix(i, colIndex) = count;
                colIndex = colIndex + 1;
                count = 0;
            end
        end
    end
    
    % If the last element in the row was 1, store the count
    if count > 0
        outputMatrix(i, colIndex) = count;
    end
end

t1count{1,k} = outputMatrix; 

end

%make 0s Nans to average
%okay now take average across rows to find avg t1 dwell time per cell
t1countavg = {};

for i=1:size(sims,1)
    in = t1count{1,i};
    in(in==0)=NaN; 
    mnz = nanmean(in');
    t1countavg{1,i} = mnz; 
end


%average of shape indices per group over each timepoint 
siavg = {};

for j = 1:size(sims,1) %groups of sims (passive, active)

            use = times{1,j};

            avgs = nanmean(use);
            siavg{1,j} = squeeze(avgs);
end

